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ABSTRACT 

■ Within a fully relativistic framework, we derive and solve numerically the perturba- 
tion equations of relativistic stars, including the stresses produced by a non-vanishing 

\ shear viscosity in the stress-energy tensor. With this approach, the real and imaginary 

^p H , parts of the frequency of the modes are consistently obtained. We find that, approach- 

2 ' ing the inviscid limit from the finite viscosity case, the continuous spectrum is reg- 

ularized and we can calculate the quasi-normal modes for stellar models that do not 
admit solutions at first order in perturbation theory when the coupling between the 

• i— i . 

, polar and axial perturbations is neglected. The viscous damping time is found to 

■ 

■ agree within factor 2 with the usual estimate obtained by using the eigenfunctions 

of the inviscid limit and some approximation for the energy dissipation integrals. 
We find that the frequencies and viscous damping times for relativistic r— modes lie 
between the Newtonian and Cowling results. We compare the results obtained with 
homogeneous, polytropic and realistic equations of state and find that the frequen- 
cies depend only on the rotation rate and on the compactness parameter Mj R, being 
almost independent of the equation of state. Our numerical results for realistic neu- 
tron stars give viscous damping times with the same dependence on mass and radius 
as previously estimated, but systematically larger of about 60%. 

Key words: stars:oscillations - stars:neutron - gravitational waves. 



1 INTRODUCTION 

Since the discovery of the r— mode instability (Andersson, 1998), the study of their astrophysical 
relevance has become a continuously growing field. In particular, the gravitational radiation from 
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the r— mode instability has been proposed to be the reason why all observed young neutron stars 
have rotational frequency much smaller than the Keplerian frequency (Andersson et ah, 1999; 
Lindblom at ah, 1998), to provide a way to identify strange stars as persistent sources of grav- 
itational waves (Andersson et ah, 2002), or to be the reason why accreting neutron stars in low 
mass X-ray binaries rotate in a narrow range of frequencies (Bildsten, 1998; Wagoner, 2002). A 
comprehensive review on gravitational waves from instabilities in relativistic stars, including the 
r— mode instability, can be found in Andersson (2003). 

Several issues concerning the r— modes have been object of controversy, among which we 
must mention the existence of the continuous spectrum, and the way in which viscosity can sta- 
bilize the stars. It was soon pointed out (Kojima, 1998; Ruoff & Kokkotas, 2001) that, working 
at first order in perturbation theory and neglecting the coupling between different harmonics and 
between the polar and axial perturbations, leads to the existence of a continuous spectrum and 
makes doubtful the existence of the modes. As shown later, going beyond the low frequency limit 
and the Cowling approximation does not remove the continuous spectrum (Ruoff & Kokkotas, 
2002). However, hydrodynamical numerical simulations (Gressman et ah, 2002), or evolutionary 
descriptions taking into account the polar/axial coupling (Villain & Bonazzola, 2002), seem to 
indicate that r— modes do exist, so that the continuous spectrum may be interpreted as an artifact 
due to an inconsistency of the perturbative expansion when a ~ fl With this motivation, more 
sophisticated methods to solve the eigenvalue problem in a general inviscid case have been devel- 
oped (Lockitch et ah, 2001; Lockitch et ah, 2003; Lockitch et ah, 2004). On the other hand, in any 
realistic situation viscosity is actually present. An interesting question to address is then if the in- 
clusion of viscosity in the game does affect the existence of the continuous spectrum, as suggested 
for example in Lockitch et al. (2004). It is believed that bulk/shear viscosity limit the instability at 
high/low temperatures, respectively. But the situation is complicated because the effects of super- 
fluidity in the inner core, hyperon viscosity, or the core-crust shear layer are uncertain (Andersson 
et ah, 2004; Glampedakis & Andersson, 2004; Lindblom & Owen, 2002). 

To our knowledge, in previous works the effect of viscosity in damping the r— modes has 
been studied in the following way. The star is assumed to be composed of a perfect fluid and the 
eigenvalues and eigenfunctions of the modes are computed by solving the adiabatic, perturbed 
Newtonian or Einstein's equations. The viscous damping time is then computed by evaluating 
some suitably defined integrals that express the energy content of the modes. This approach has 
been used, for example, in the seminal work about viscosity on neutron star oscillations by Cutler 
& Lindblom, (1987), and it is justified by the fact that viscosity is indeed very small. However, 
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even a very small amount of viscosity may be crucial to change the role of the r-mode instability 
in astrophysical processes, therefore this limit deserves to be studied more accurately. In this work, 
we introduce the shear viscosity in the stress-energy tensor from the beginning, and we evaluate 
consistently both the frequency and the damping time of the mode. By this approach, one could 
determine more precisely the window of instability, for a given stellar model, provided that the 
microphysical input (EOS, viscosity) is known. A similar self-consistent inclusion of the heat 
transfer corrections to the polar perturbations has recently been done (Gualtieri et ah, 2004). An 
interesting results of this work is that shear viscosity, even if very small, regularizes the continuous 
spectrum. 

The paper is organized as follows. In Sec. 2 we derive the equations of stellar perturbation 
with shear viscosity. In Sec. 3 we discuss the numerical integration of the perturbed equations 
and the boundary conditions. In Sec. 4 we examine the results for the different stellar models 
(homogeneous, polytropic, realistic) and in Sec. 5 we draw the conclusions and discuss the domain 
of applicability of our approach and future extensions. In the Appendix, an analytic solution for 
Newtonian r— modes with viscosity is derived. 



2 THE PERTURBED EQUATIONS WITH SHEAR VISCOSITY 

We consider a star rotating uniformly with angular velocity VL. At first order in Vt (or, more pre- 
cisely, at first order in the rotational parameter e = Vl/VL N with Vt N = ^2M/R 3 ), the stationary 
background is described by the metric (Hartle, 1967; Hartle & Thorne, 1968) 

ds 2 = g^dx< x dx v = -e u{r) dt 2 + e x(r) dr 2 + r 2 (d$ 2 + sin 2 $d V 2 ) - 2r 2 u{r) sin 2 0d(pdt , (1) 

where uo(r) represents the dragging of the inertial frames. It corresponds to the angular velocity of 
a local ZAMO (zero angular momentum observer), with respect to an observer at rest at infinity. 
The 4-velocity of the fluid is simply = (e~ u / 2 , 0, 0, VLe~ v l 2 ), and the stress-energy tensor is 

i2» = (p+p)«M 0) + w s». (2) 

We assume that viscosity does not affect the stationary axisymmetric background, because the 
shear tensor and the expansion do vanish there. Therefore, in this regime the Einstein field equa- 
tions reduce to the standard TOV (Tolman-Oppenheimer-Volkoff) equations plus a supplementary 
equation for the frame dragging uo(r): 

uj )TT - ^4vr(p + p)e x r - ^ uj )T - 16vr(p + p)e x u = (3) 
where u>(r) = Q —u(r). 
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2.1 Perturbations with axial parity 

In this paper we focus on purely axial perturbations of the metric ([T]). The perturbed metric, ex- 
panded in spherical harmonics and in the Regge- Wheeler gauge, depends only on the axial vector 
harmonics 

t=(-^sin^) (4) 

where a = $,(p and Y 1 ™^, ip) are the scalar spherical harmonics (in the following we omit the 
harmonic indexes lm). The Einstein equations can be shown to reduce to a set of equations for two 
metric variables (h , h\) and for the axial component of the fluid velocity (Kojima, 1992), whose 
harmonic expansions read: 

hoa(t,r,#,<p) = ho{r,a)S a {$,ip)e- i<rt 
h la (t,r,#,v) = hiMSa&ipy-™ 

6u a (t,r,$,v) = ^Z(r,a)r b S b ($,v)e- iat . (5) 

Here, 7^ =diag(l, sin 2 $) is the metric on the two-sphere, a is a (complex) frequency, 5u r = 0, 
and 5u l can be obtained imposing u^u^ = — 1. The axial perturbations are then described by only 
three degrees of freedom: two, ho, hi, for the metric, and one, Z, for the fluid. Notice that for the 
metric variables we follow the convention used in Kojima (1998) and Ruoff & Kokkotas (2002), 
but we use a different fluid variable. Our variable Z is related to the u 3 in Ruoff & Kokkotas (2002) 
(indicated with RK) by 

Z = e^' 2 u* K - h (6) 

and our functions A, v are related to those used in Ruoff & Kokkotas (2002) by A = 2\ RK , v = 
2u RK . In order to simplify the form of Einstein's equations we define the following variables 

K = e (x - u)/2 h (7) 
V = -ie (l/ - A)/2 /n (8) 



H = e 



-A 



K' + ( - — - - - ) K — (a — muj)e x - v V 



(9) 



2 r 

which are related to those of Ruoff & Kokkotas, (2002) by 

K = K RK 
V = -iV* K 

H = K RK . (10) 
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Notice that we have introduced an auxiliary variable, H , in order to have, at least in the perfect 
fluid case, a system of first order differential equations. 



2.2 Einstein's equations 

Taking into account viscosity, the stress-energy tensor has the form (Misner et al, 1973) 



r vv = (p + p)u^u„ + pg^ - 2-qa^ - 2(g lxu u a . a . 



(ID 



Here 77, ( are the shear and bulk viscosity coefficients (see Cutler & Lindblom, (1987) for a dis- 
cussion on the meaning of these coefficients), and 



with P^ being the projector onto the subspace orthogonal to 



(12) 



(13) 



In the interior of the star, after expanding the perturbed tensors in spherical harmonics, Einstein's 
equations perturbed at first order are: 

V-A' 2" 



K' 
V 



K + (a- mio)e x ~ u V + e x H 



(14) 



mr 

T 

m 



-(a - muj)K + ^- [uj'H - 16tt( P + p)u (K + Ze^ )/2 )} 
+16nir] 



e ^Z - 'j-aur 2 (Ze x ^ + Ke^ 2 



(15) 



H' = 



In + ^-^K + 16n(p + p) [Ze^ 2 + K] - 2^V W 



- 167Tir]mQ^j^e x/2 - u Z 



(16) 



-(a — mu)H 
X 



_ 2mu\ „ 
(T-msH — i Z 



= — — ^ + — — e~ A AT - 167ri?7X 
r 2 A 

= e" A / 2 (z' - 2 -Z^ + (a - m^e^/V 

p {v-X)/2 

= -{a - mQ)e (u - x)/2 K + if — x 



(17) 
(18) 



(p + p)r 2 



x 



(r 2 r)X)' - r)(A - 2) e x ' 2 Z - (e^Z + e^ 2 K) 



A 



(19) 
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where A = 1(1 + 1) and we have introduced the new variable X (Eq. (fTSl ). We note that this quan- 
tity comes directly from the expansion in axial vector harmonics of the shear tensor components 
a r $, a rip . Having written the equations, we would like to make some relevant remarks. 

• We have neglected the coupling between axial perturbations with harmonic index I and polar 
perturbations with harmonic indexes I ± 1, as in Ruoff & Kokkotas (2002) and Kojima (1993). 
The reason is that in this paper we focus on the effect of viscosity on r-mode oscillations, and we 
want to separate different effects. We point out that the coupling terms between polar and axial 
perturbations affects the shape of the continuous spectrum (Ruoff et ah, 2002). However, in this 
paper we show that, having neglected this coupling, viscosity removes the continuous spectrum. 

• Under the approximation described above, bulk viscosity is not coupled to axial perturbations, 
because it enters into the equations only through the axial-polar I I ± 1 couplings. Thus, we will 
only study the effects of shear viscosity, leaving the investigation of the effects of bulk viscosity 
for future work. 

• In the perfect fluid case r] = 0, we have three first order differential equations, (fl4l4T6T) . and 
three algebraic relations, (fl7l - fT9l . Notice that H can be given, alternatively, by the differential 
equation (fToT) or by the algebraic relation (fTTt . While the inviscid limit of equations dT~4l - fT31) and 
(fl7l - fT9l has been derived before (e.g. Ruoff & Kokkotas, (2002)), the differential equation for H 
(fToT) has never been derived so far; we have checked numerically that it is equivalent to determine 
H by (fTTt or by integration of (fToT) . 

• In the perfect fluid case r] = the equation (fT9l) for Z is a simple algebraic relation, and Z 
can be found by dividing by (cr — mVt + ^jp) ; this is the origin of the continuous spectrum which 
appears when this term vanishes. When viscosity is present, extra terms do appear involving the 
second derivative of Z so that, in order to solve our equations in the frequency domain, we prefer 
to recast Eqs. (fT8l - fT9l as 

Z' = -Z + e x l 2 X-(a-mti)e^l 2 V (20) 
r 



A v ; 



{ {P+P) c (X-vV2 



L-mn+ ^) Z + (a- mQ)e^ 2 K 



. , (2D 

V 

Because of numerical problems associated to the 1/rj term, we cannot integrate this system of 
equations for r] < 10~ 6 km -1 . 1 



1 Usually the typical viscosity at T = 10 7 K and p = 10 15 g/cm 3 is r\ ~ 10 23 g/cm/s = 2.4 • 10 11 km x . 



Relativistic r -modes and Shear viscosity: regularizing the continuous spectrum. 1 

In the exterior of the star, the stress-energy tensor vanishes and the equations for the metric 
perturbations are 

K' = ( A' + - \k + (a - muj)e 2X V + e x H (22) 

,2 



V 

H' 
H 



— (a — muo)K + 



mr 



A 



r r z A 

A 2 2muo' 

— ^~ V + — ^ e K 

a — muj r z A 



(23) 
(24) 
(25) 



where we can alternatively use (l24"|) or d25t to determine H. 



2.3 Cowling approximation and Newtonian limit 

Since in Section|4]we shall make a detailed comparison with previous works, it is useful to write 
the equations in the Cowling approximation and in the Newtonian limit. The equations in the 
relativistic Cowling approximation can easily be obtained. By neglecting all metric perturbations, 
and considering as dynamical equations of the system ST^ = 0, we end up with a single second 
order differential equation for the fluid variable Z, which we write as a system of two first order 



equations: 

21 = - 
r 

X' = - 



Z + e x/2 X 



(26) 



2 rj 
- + — 

T Tj 



x + 



(A -2) 



1 - 



r raauoe 



A 



V V Ay 

Notice that in the limit 77 — > 0, it reduces to 

2mu 



a — mVt 



A 







(27) 



(28) 



showing the well known problem of the continuous spectrum due to the dragging of the inertial 
frames. The Newtonian limit can be obtained by taking e x = 1, e u = 1, pushing the speed of light 



to infinity (so that in our units p/p 
resulting equations are: 



0) and neglecting the relativistic corrections (cu = fi). The 



Z' 
X' 



-z + x 

r rj 



r 2 maQ 


z-.e( 




I 1 A 


77 V 



A-2 



a 



A 



-mQ Z 



(29) 
(30) 



with the very well known inviscid limit o 



(A-2) 



As shown in the Appendix, equations (1291430b admit an analytical solution for models with 
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constant p and 77. Furthermore, with a suitable change of variable, they can be integrated numeri- 
cally for any value of r\. 



3 NUMERICAL IMPLEMENTATION 

The system of equations (fl?14T6T) . d27)l ). (OTT) coupled to the TOV equations and Hartle's equation 
© is solved by means of an adaptive step, third order Runge-Kutta integration scheme. After 
integrating the equations in the interior of the star, we integrate the exterior equations d22l - l25l) 
backward starting with an expansion in terms of the purely outgoing solutions at a large radius, as 
explained in Sec. 13.31 By evaluating the Wronskian of the interior and exterior solutions, we can 
obtain the quasi-normal mode frequency by monitoring when the modulus of the Wronskian van- 
ishes. We have implemented a Newton-Raphson scheme to find the correct frequency at which the 
Wronskian vanishes starting with an initial guess. In the inviscid case, there is only one indepen- 
dent solution that is regular at the origin; in the viscous case there are two independent solutions 
satisfying the boundary conditions at the center, as we show below. Thus, an additional boundary 
condition will be required at the star's surface. 



3.1 Boundary conditions at the center 

Expanding equations (fl4l - fT9t near r = 0, we find the following behaviour: 

K = r l+1 

V = Ar l+2 

H = (l-l)r l 

Z = Cr l+1 + Dr l+3 

where 



(31) 



A -- 

C : 

and 

A : 

B ■ 



1 + 2 



[—{a — mu) c ) + 16nir]C] 



-(a - rafi)e" c / 2 + 



in 



Pc+p, 



■A 



(a — mfl 



2rau) c 



) + 



in 



-B 



' 4 (a - mQ)(o - moj c ) + maO c —-^ + 2e Uc/2 (2l + 3)D 



1 + 2 

Uc/2 



A 



(/ - 1)(2Z + 5)^7T Pc - 167Ti77^|^" 2 <vr- 



( "(a — mQ) — 2mauj c e 



A-2~ 
A 



(32) 
(33) 
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This last relation couples the constant D to the others. The reason why we need the term propor- 
tional to r /+3 in the expansion of Z is that, keeping only the leading term in r z+1 , Eq. (fl"9l results 
in a trivial identity that does not determine the arbitrary constants. Notice that when rj = we 
have simply, 

C = - i-~^y c/2 (34 ) 

and D is no longer coupled to the other constants, so that, for a given a, there is only one solu- 
tion of the equations satisfying the conditions at the center. When r/ ^ 0, instead, there are two 
independent solutions which satisfy the conditions (|3"T1) at the center. 

3.2 Boundary conditions at the surface 

As discussed above, in order to integrate the perturbed equations for r] ^ 0, an additional condition 
needs to be imposed at the surface of the star. We impose the Israel matching conditions (Israel, 
1966), that is, we require continuity of the induced metric and of the extrinsic curvature on the 
timelike surface which separates the interior from the exterior of the star. From the continuity 
condition of the induced metric we learn that h is continuous across the stellar surface 

[h ] = 0, (35) 

where we use the notation [/] = lim e ^ (f(R s + e) — f(R s — e)) . Imposing the continuity of the 
extrinsic curvature we find that 

N = 0, [h' ] = 0. (36) 

These conditions are analogous to those found in Appendix B of Price & Thorne, (1969) for the 
polar case. When applied to our variables, they translate in 

[K]=0, [H] = 0, [V] = 0. (37) 

Requiring that the interior equations (TB14T91 match with the exterior equations (I22TE51 ) leads to 
the stress-free condition r]X = at the stellar surface, which arises directly from equations (flTl) . 

The physical meaning of the stress-free condition can be understood by looking at the New- 
tonian case (or at the Cowling approximation), in which case one cannot rely on conditions that 
involve metric variables. Consider the stress acting on an element of fluid inside a spherical thin 
layer bounded by the the surface of the star and a sphere of radius R s — e, and by two lateral 
surfaces defined by -&i =const. and $ 2 =const. The stress on the upper (r = R s ) boundary is zero 
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since there is no fluid outside. If e — > 0, the stress through the lateral surface goes to zero as e and 
the only stress is through the interior surface. But, since the volume of the fluid goes to zero with 
e, the stress acting on the interior surface must vanish, otherwise the difference would result in an 
infinite acceleration of the thin layer. In the background, the only non vanishing radial component 
of the stress tensor <Sjj (i.e. the pressure) vanishes at the stellar surface; the same condition must be 
imposed on the perturbed stress. Assuming that the perturbations of radial velocity and pressure 
are zero (purely axial perturbations), the radial components of the perturbation of the stress tensor 
in spherical coordinates are 

d ( Sv#\ rr , d ( Sv v 



SS rr = , 5S r # = rjr— , 5S rip = r/r— , (38) 

or \ r J or \ r 

and imposing their vanishing at r = R s , we find 

d ( 5v#\ n d ( 8v v 



dr\rj ]r=Ra dr \ r J ]r=Rs 

After expanding in spherical harmonics and using eq. ©, these conditions are equivalent to 



\ r 



Z = rjX = . (40) 

This equation shows that, at the boundary between two fluids, normal shear stresses must be con- 
tinuous. 

In order to better understand the relationship between the Newtonian condition d40l and the 
relativistic condition r)X(R s ) = 0, we note that the shear tensor components cr ri? , a rip , expanded in 
axial vector harmonics, give exactly the perturbation function X, defined in (fT8t : so our matching 
condition T]X(R S ) = is equivalent to the vanishing of r/a r ^, rja rip at the stellar surface. 

On the other hand we recall that the Israel matching conditions, as shown in Israel (1966), are 
equivalent to impose that 

[Gra]=0 O=0,. ..,3) (41) 

which, from Einstein's equations, are equivalent to T ra (R s ) = 0. The a = a = components 
of this equation give 

T ra (Rs) = ((P + P)u r Ua + P9ra + Wra) r= R s = {Wra)r=R s = > ( 42 ) 

where we have used the fact that, for a purely axial perturbation, u r = + 8u r = and pg ra = 
on the stellar surface. 
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3.3 Outgoing wave condition at infinity 

After imposing the condition r/X = at r = R s , for each value of the (complex) frequency a the 
perturbed equations admit only one solution extending up to radial infinity. In practice, we cannot 
extend our integration up to a large distance from the center because the amplitude of the outgoing 
component is exponentially increasing 2 . However, we found it sufficient to integrate backward, 
starting from r = 200 R s with an asymptotic expansion of the outgoing wave up to first order in 
1/r: 



rout _ :~fi\fi\ „i&r* i n I ^ 



yout = ir j 1 + J e «n-. +0 | _ ] (43) 

with 



a 



3kr - 2 
icr(l + icr) 



icr(l + icr) 

Notice that V are proportional to r because h ,hi, in the Regge-Wheeler gauge, are propor- 
tional to r; however this apparent divergence is not a real problem, being only an artifact of the 
gauge choice (see for example Thorne & Campolattaro, 1967). 

The Wronskian W of the two solutions (interior and exterior) of the system of equations is 

r — 2M 

W = — r (KV out - VK out ) . (45) 

H(<7 — muj) 

We have numerically checked that W is constant independently of the radius at which it is evalu- 
ated (typically with 7-8 significant digits). 



4 RESULTS 

We shall now discuss the results of the numerical integration of our equations. For comparison 
with previous works and for testing purposes, we begin discussing the case of constant density 
stars, later showing the results for polytropic stars and realistic neutron stars. In each case, we 
compare the Newtonian limit, the Cowling approximation, and the relativistic calculation, with 
the purpose of establishing the qualitative differences between the three approaches and to check 
whether our results converge to known results in the inviscid limit. 



2 See Kokkotas & Schmidt, (1999) for a discussion. This problem can be overcome by extending r to the complex plane (Andersson et at, 1995), 
or by deriving a solution for the equations outside the star in terms of recurrence relations (Leaver, 1985; Leins et at, 1993). 
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4.1 Homogeneous Stars 

For constant density stars the frequencies of the r— modes have been shown to lay outside the 
continuous spectrum (see e.g. Kokkotas & Stergioulas, 1999 for the full GR results). We consider 
a uniform density star with a central energy density of 10 15 g/cm 3 , mass of 1.086 M Q and radius 
of 8.02 km (M/R = 0.2). We choose the rotational parameter e = 0.3. In Fig. 1 we show the 
r-mode frequency versus the viscosity parameter rj, assumed to be constant throughout the star. 
The Newtonian, Cowling, and General Relativistic results are shown by dotted, dashed, and solid 
lines, respectively. The Newtonian and GR frequencies in the inviscid limit (1064 Hz and 1144 
Hz) are indicated as dotted and solid horizontal lines, while the shadowed region indicates the 
continuous spectrum. Fig. 1 clearly shows that as the viscosity decreases the inviscid Newtonian 
and GR results are recovered, while (as expected) the frequency in the Cowling approximation 
falls inside the continuous spectrum. Since the convergence to the inviscid limit is reached for 
rj ks 10~ 5 km -1 , the mode frequency in the Cowling approximation can be found by extrapolating 
the dashed line for rj — > 0; the corresponding value is 1244 Hz, about a ten percent larger than the 
GR value. It should be stressed that the mode frequency in the Cowling approximation had never 
been obtained before for perfect fluids. 

Since for numerical reasons we must restrict our calculation to values of rj > 10~ 6 km -1 , 
we are always in the region where the viscous damping dominates the instability. Therefore, the 
inverse imaginary part of the frequency which we show in Fig. 2 is a damping, not a growth 
time. Two interesting features are shown in Fig. 2. Firstly, the Newtonian, Cowling and GR damp- 
ing times differ for less than 20% because the metric corrections are small. Notice that in the 
Newtonian approach the stellar background is obtained solving the relativistic equations and the 
perturbations are studied using the Newtonian equations of hydrodynamics. Secondly, all damp- 
ing times show the expected l/rj behaviour. More precisely, we can compare our results with the 
analytic back-of-the-envelope formula of Cutler & Lindblom, (1987) for the dissipative time scale 
of the shear viscosity, based on a quasi-uniform density Newtonian model: 



T (/-1)(2/ + 1)?/ ( } 

For our model r = 3 x lO -8 /^, with r\ in km -1 and r in s. This estimate is also shown in Fig. 
2 (thin solid line) and it is, surprisingly, in better agreement with the GR results than with the 
Newtonian ones. 
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It is useful to introduce the Ekman number, which can be defined as 



where P is the period of the mode. This number represents the ratio between the viscous term 
and the Coriolis force. Typically, it is smaller than 1CT 7 (Cutler & Lindblom 1987). For numerical 
reasons, our calculations cannot be extended to values of E s smaller than 0.1, which is still far 
from real models. Nevertheless, it allows to obtain the r— mode frequencies in the inviscid limit 
with a good accuracy. 

4.2 Polytropic Stars 

We now turn to polytropic stars. For these models, and within the same approximations we use to 
derive our equations, the r— mode was found to disappear for certain ranges of the polytropic index 
n (Ruoff & Kokkotas, 2001; Ruoff & Kokkotas, 2002). Only for very compact stars n < 0.8, the 
r— mode frequency lies outside the continuous spectrum and could be found. We have considered 
a polytropic model with n — 1, and the same compactness and rotation parameter as in previous 
section (M/R = 0.2, e = 0.3). This model has mass M = 1.74M© and R = 12.86 km. As shown 
in Ruoff & Kokkotas, 2001, in the inviscid case we do not find the r— mode because it lies inside 
the continuous spectrum. In Fig.|3]we show the results obtained when viscosity is included in the 
calculation. In the Newtonian case (dotted line), the inviscid limit is nicely recovered as before. As 
for the Cowling (dashed line) and RG (solid line) calculations, we can follow the r— mode inside 
the continuous spectrum until convergence to the inviscid limit is reached. The corresponding 
damping times are shown in Fig.El Similarly to constant density models, the Relativistic damping 
time lies between the Newtonian and the Cowling calculation, and they agree within a factor 2. 
For comparison we also include the estimate given by Eq. d46l) using the average density of the 
star, which overestimates the damping time of about 50 %. Notice that using the average density in 
Eq. d4oT) . the damping time (for models with constant rf) depends only on pR 2 oc M/R, therefore 
it gives the same result for stars with the same compactness. 

The results shown in this section can also be compared with a recent work by Villain & Bonaz- 
zola, (2002), in which the coupling of the polar and axial perturbations was fully taken into ac- 
count, and the perturbed Euler equations were integrated in the Cowling approximation with spec- 
tral methods. For a similar polytropic model they found an r— mode with a frequency about a 10% 
higher than the Newtonian value, in good agreement with our results. 

To complete our discussion about the existence of r-modes in the continuous spectrum, in 




(47) 
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Fig. |5] we show the behaviour of the real part of the r— mode frequency (er) as a function of the 
polytropic index n for models with compactness Mj R = 0.2, period of 1 ms, and a shear viscosity 
rj = 10~ 6 km -1 . The period and compactness have been chosen to allow for direct comparison 
with the results of Ruoff & Kokkotas (2001) who could not find the r-modes for n larger than a 
certain value, when the real part of the frequency reached the continuous spectrum. By introducing 
a small amount of viscosity, the frequency can be calculated for all polytropic indexes, or for any 
other stellar model, even if we stay at the most basic level of approximation: first order in the 
rotation parameter and neglecting the coupling between the axial and polar parts. Note that when 
the mode lays outside the continuous spectrum we obtain results very similar to Ruoff & Kokkotas 



Some attention must be paid to the fact that viscosity breaks the degeneracy of the r— modes 
and allows for the existence of a family of modes with an arbitrary number of nodes. A simple 
analytical solution for homogeneous stars in the Newtonian case and for constant 77 is detailed in 
the Appendix. In this particular case, the eigenfunctions are simply given in terms of the spherical 
Bessel functions Z = r ji(kr). For the / = 2 case we have 



Indeed, we can obtain a family of solutions labeled by different wave-numbers k. In Fig.|6]we show 
the first three eigenfunctions of the Newtonian r modes for constant density stars and for polytropic 
stars plotted with dashed and solid lines, respectively. The numerical results for homogeneous stars 
have been checked to coincide with the above analytical solution. All eigenfunctions have been 
normalized to their surface value. The analytical solution shows that, for low viscosity, the leading 
order correction to the frequency is proportional to a (r]/p) 2 — iri/p, with cr = ^^mVl. The linear 
dependence of the inverse damping time with r/ has already been discussed in the previous section; 
here we notice that the correction to the frequency has a quadratic dependence on 77. Remarkably, 
we find that even for relativistic perturbations, the numerical results can be well fit by this quadratic 
behaviour with a very good accuracy. 

A last relevant result is that in the fully relativistic calculation we could not find more solutions 
than the fundamental mode. The other overtones found to exist in the Newtonian or Cowling 
approaches, seem to be incompatible with the pure outgoing wave condition imposed at infinity. 
In Fig. [7] we show a comparison of relativistic r— mode eigenfunctions for constant density stars 
(for which the inviscid limit is available) with different viscosities: 77 — (solid), 77 = 1.4 x 10~ 7 
km -1 (dots), 77 = 10~ 6 km -1 (dashes), and 77 = 10~ 4 km -1 (dash-dot). The figure shows that the 



(2001). 




(48) 
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eigenfunction of the fundamental r— mode approaches continuously the inviscid limit, which is 
close to the r 3 dependence. It should be reminded that both in the Newtonian and in the Cowling 
approach the gravitational perturbations are neglected; thus, the fact that the relativistic r— modes 
do not admit solutions with large wavenumbers may be an indication that they only exist for fluids 
where the axial fluid motion is totally decoupled from gravitational perturbations. 



4.3 Realistic Neutron Stars 

In this section we consider stellar models constructed with realistic EOSs and realistic viscosity 
profiles with the purpose of establishing if any imprint of the equation of state could be carried 
by the r— mode gravitational wave spectrum. At low density (below 10 12 g/cm 3 ) we use the BPS 
(Baym et ah, 1971) equation of state, while for the inner crust 10 12 < p < 10 14 g/cm 3 we employ 
the SLy4 EOS (Chabanat et al, 1998). At high density (p > 10 14 g/cm 3 ) we will consider two 
different EOSs of neutron star matter representative of the two different approaches commonly 
found in the literature: potential models and relativistic field theoretical models. As a potential 
model, we have chosen the EOS of Akmal-Pandharipande-Ravenhall (Akmal et ah, 1998, hereafter 
APR). As an example of the mean field solution to a relativistic Walecka-type Lagrangian we 
have used the parametrization usually known as GM3 (Glendenning & Moszkowski, 1991). For a 
detailed discussion and comparison of many equations of state from the two families see the recent 
review by Steiner et al. (2005). Setting the compactness parameter to M/R = 0.2, APR gives a 
mass of 1.53 M and a radius of 11.3 km, while GM3 gives M = 1.72M and R = 12.7 km. 
For a polytropic EOS with n = 1 the corresponding mass and radius for the same compactness 
are M = 1.74 M Q and R = 12.8 km. For cold neutron stars below 10 9 K, neutrons in the inner 
core become superfluid and the dominant contribution to the shear viscosity is electron-electron 
scattering (see e.g. the review by Andersson & Kokkotas, 2001); in this regime 77 can be written as 

p ee = 6 x 10 18 p?5 T 9 _2 g/cm/s = 1.48 x 10" 15 p? 5 T 9 " 2 km" 1 (49) 

with p 15 and T 9 being the density and temperature in units of 10 15 g/cm 3 and 10 9 K. Having this 
in mind, in this section we will use a viscosity coefficient with a quadratic dependence on density 

77 = 7/0^ km" 1 , (50) 

where p c is the central density. Since old neutrons stars are nearly isothermal, we will parametrize 
our results as a function of the constant 770 that includes the temperature dependence. 

In Fig. [8] we show the r— mode frequency as a function of 770, comparing the three EOSs: 
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polytrope (dots), APR (dashes) and GM3 (solid lines). In all cases the rotation period is P = 2 
ms. As we can see in the figure, the frequency is rather insensitive to the particular details of the 
EOS, provided that the rotation frequency and M/R are the same. The Newtonian limit depends 
only on the angular frequency (0) and the relativistic correction enters through the frame dragging. 
Since this correction goes as u;/f2 « I / R? oc M/R, the leading order contribution to the frequency 
is a function only of the compactness. 

In Fig. |5] we show the damping times for the same models as in Fig. [5] For constant density 
stars, the « pR 2 dependence of the damping time translates into a Mj R dependence for models 
with constant r\. In realistic, cold (T < 10 9 K) neutron stars, the density is not constant and the 
viscosity is dominated by the electron-electron scattering process (l49t . Thus, the analytic result 
(146b underestimates the viscous damping time. An improved calculation for n — 1 polytropes 
taking into account the density profiles and using the shear viscosity (l49t (Andersson & Kokkotas 
2001) gives 
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2.2 x 10' ( — w ) ) T*s. (51) 
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Assuming a general dependence of the viscosity of the form of Eq. (l50l . and using Eq. d49hhe 
above damping time can be shown to satisfy 

= (UMe\ (R\ 5 / Pc \ /3.26X10-WN 
1 M ) VlOkmy/ Vl0 15 g/cm 3 ; V Vo J 



Our relativistic calculations show approximately the same dependence on mass and radius as Eq. 
d52l . but the damping times are systematically larger of about 60%. We found that a better fit for 
the realistic neutron stars (APR, GM3) as well as for the n — 1 polytrope is given by 

/L4Mq\ / R \ s ( Pc \ / 5.22x lO^km- 1 ^ g 

V M J VlOkmyl Vl0 15 g/cmV V % / S ' 

In Fig.|9lwe show, together with the numerical results (thick lines), the results corresponding to 
the previous fit (thin lines). The good agreement between them is apparent. 



5 FINAL REMARKS AND CONCLUSIONS 

We have investigated the r— mode spectrum within a perturbative relativistic framework consis- 
tently including the effects of shear viscosity. A first interesting result of our study is that the 
continuous spectrum, which was an artifact of the level of approximation (first order in rotation, 
neglecting completely the coupling between polar and axial perturbations), can be regularized and 
the frequencies and viscous damping times of the modes can be calculated. This can be understood 
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if one analyzes the real cause of the existence of the continuous spectrum. By considering only 
purely axial perturbations, and working at first order in the rotation parameter, each spherical fluid 
layer is effectively decoupled from the neighbour (indeed, the perturbations of density, pressure, 
and radial velocity are zero). Therefore, each layer can oscillate independently and have a differ- 
ent characteristic frequency. This effect can be produced by differential rotation or, in a relativistic 
framework, by the dragging of the inertial frames. It has been discussed in the literature that in- 
cluding the coupling up to some level, the shape of the continuous spectrum changes (Ruoff & 
Kokkotas, 2002), but still one cannot always find a mode. In principle, by fully including the cou- 
pling (see e.g. Lockitch et ah, 2001; Lockitch et ah, 2003; Villain & Bonazzola, 2002), one can 
look for a global mode that oscillates with a unique characteristic frequency, since different layers 
feel each other and their motion is coupled. We have produced the same effect by introducing 
shear viscosity. The stress between the different layers couples them and results in the existence of 
a unique global r— mode, instead of a continuous spectrum. Therefore, even working at the most 
crude level of approximation (purely axial perturbations, first order in rotation) one can calculate 
the r— modes (complex) frequency for every stellar model. 

Another interesting outcome of our study is that, in the regime where viscosity dominates, the 
damping times of relativistic modes are about 40-60% larger than those obtained in the Newtonian 
case. For homogeneous stars the differences are smaller, and the semi-analytic estimate obtained 
by Cutler & Lindblom (1987) fits well our numerical results. For realistic neutron stars we find 
that our relativistic results give damping times systematically larger (about a 60%) than previous 
Newtonian estimates (Andersson & Kokkotas, 2001), but with the same dependence on mass and 
radius. An accurate fit to realistic neutron stars, valid also for n = 1 polytropes gives 



other in such a way that the previous estimate is actually closer to the results for Newtonian 
constant density stars (Kokkotas & Stergioulas, 1999; Andersson & Kokkotas, 2001) than the 
more elaborated calculation for n — 1 polytropes (Andersson & Kokkotas, 2001). 

Therefore, we have shown that by including the shear viscosity in the stress-energy tensor we 
can obtain an accurate estimate of the damping time of the r— mode in neutron stars, provided that 
realistic equation of state and viscosity are known. We were able to extend our results to arbitrarily 
low viscosity in the Newtonian case, but further improvements, such us including higher order 
couplings or devising better numerical approaches to deal with smaller values of 77, are needed 




(54) 



Notice that the corrections due to general relativity and the internal structure almost cancel each 
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before we can compute accurately the window of instability of the r— modes in a fully relativistic 
framework. Notice that axial-polar coupling (second order in rotation) must also be included to 
take into account the effect of bulk viscosity. 

In this paper we have not considered the damping associated to the viscous boundary layer 
in the core-crust interface, which is thought to be one of the most efficient mechanisms to damp 
the r— mode oscillations (Bildsten & Ushomirski 2000; Rieutord 2001). If the crust is assumed to 
be rigid, the problem can be formulated as an Ekman problem in a spherical layer (Glampedakis 
& Andersson, 2004). In practice this means to substitute the boundary conditions at the stellar 
surface by a "no-slip" condition at the core-crust interface, i.e., Z(r = R c ) = 0, and to include the 
corrections due to the perturbation of the pressure. Shear viscosity in this case is essential since it 
ensures that the perturbations of the velocity deviate from the inviscid solution and go to zero in 
a thin layer near the crust (the Ekman layer) with characteristic width ^ m \^W S . Recent results 
(Glampedakis & Andersson, 2004) show that the damping rate and the correction to the frequency 
due to the presence of the Ekman layer go as (Q/ ' R c ) l l 2 , with a numerical factor of the order of 
unity. This issue will be addressed in a future work. 
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Figure 1. The r— mode frequency is plotted as a function of the shear viscosity parameter rj for a uniform density model with M/R = 0.2 and 
rotational parameter e = 0.3. The Newtonian, Cowling and GR results are shown, respectively, with dotted, dashed, and solid lines. The thin 
horizontal lines indicate the corresponding Newtonian and GR inviscid limits, while the continuous spectrum is the shadowed region. 




Figure 2. The r— mode damping time is shown as a function of the shear viscosity parameter rj for the same model as Fig.Q The Newtonian, 
Cowling and GR results are shown, respectively, with dotted, dashed, and solid lines. The thin solid line refers to the simple estimate obtained using 
Eq. 
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Figure 3. The r— mode frequency is plotted as a function of r) for a polytropic star with n = 1, M/R = 0.2 and rotational parameter e = 0.3. 
The Newtonian, Cowling and GR results are shown, respectively, with dotted, dashed, and solid lines. The horizontal dotted line indicates the 
corresponding Newtonian inviscid limit, while the continuous spectrum is the shadowed region. 




Figure 4. The r— mode damping time is plotted as a function of r\ for the same polytropic star of Fig.|5] The Newtonian, Cowling and GR results 
are shown, respectively, with dotted, dashed, and solid lines. The thin solid line refers to the simple estimate obtained using Eq. EU. 
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Figure 5. The real part of the r— mode frequency is plotted (a) as a function of the polytropic index n for models with a shear viscosity parameter 
rj = 10~ 6 km -1 , compactness M/R = 0.2 and period of 1 ms. 
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Figure 6. Newtonian r— mode eigenfunctions for constant density stars (dashed lines) and polytropic stars (solid lines) for rj = 10~ 6 km -1 . The 
numerical results corresponding to constant density stars are found to be in agreement with the analytic eigenfunctions derived in the appendix. 
The eigenfunctions of polytropic stars have the same qualitative shape. 
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r/R 

Figure 7. Comparison of relativistic r— mode eigenfunctions for constant density stars and different viscosities: rj = (solid), r\ = 1.4 X 10 — 7 
km -1 (dots), ?7 = 10~ 6 km -1 (dashes), and r) = 10 -4 km -1 (dash-dot). In the relativistic case, no mode solution with nodes could be found 
since only the fundamental r— mode admits a pure outgoing wave solution. 
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Figure 8. Comparison of the r— mode frequency as a function of the the viscosity coefficient at the center, rjo, for a polytropic star with n = 1 
(dotted line), the APR (dashed line) and the GM3 (solid line) equations of state. The compactness parameter in all cases is M/R = 0.2 and the 
rotation period is 2 ms. 
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Figure 9. Comparison of the r— mode damping time as a function of r?o for a n = 1 polytrope (dotted line), the APR (dashed line) and the GM3 
(solid line) equations of state, for the same models considered in Fig.|8] The thin lines, whoch show the estimates given by Eq. are difficult to 
distinguish from the real results. 
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APPENDIX A: AN ANALYTICAL SOLUTION FOR THE NEWTONIAN i?-MODES 



The Newtonian equations (I29H30I) describing the axial perturbations can be written as a single 
second order ODE 

r 2 Z" - 1(1 + l)Z + r 2 (a a + iUo - a ^j Z = (Al) 

where a is the (real) frequency in the inviscid limit a = ^j^mVt. If we assume that v = - p is a 
constant, this is just a version of the Riccati equation, 

z 2 y" + (k 2 z 2 - 1(1 + 1)) y = (A2) 
which has solutions 

y = kz(Aj l (kz) + By l (kz)) (A3) 

where and yi(kz) are spherical Bessel functions of the first and second kinds. In our problem 

we impose regularity of the function at the origin, which implies B = because yi(kz) diverge in 
z = 0. Therefore, we are left with the simple (complex) solution 

Z = Ax 3l (x) (A4) 

where x = kz and k 2 = <j <j + -(o — cr ). Notice that this functions have the correct r l+1 be- 
haviour near the origin. The eigenmodes and eigenfunctions are then fixed by imposing boundary 
conditions at a given radius r = R. Consider the / = 2 case and let A = 1, the solution is: 

( 3 \ 3 

Z(x) = I — — 1 ) sin x cos x . (A5) 



The boundary condition Z' — 2Z /r = at r = R implies 

(5x 2 — 12) sinx — x(x 2 — 12) cosx = (A6) 

which has only real solutions. The first few zeros are approximately x = 2.5, 7.14, 10.5, 13.8, 17.0 . 
The frequency of the fundamental r— mode corresponds to kR = 2.5, but there is an infinite 
number of higher order overtones. For each wavenumber, the deviation from the inviscid limit 
(s = a — cr ) can be readily calculated 

s = ^f^{a ^-iu). (A7) 

Notice that in the low viscosity limit, the damping time goes exactly as rj/ pR 2 , with a multiplica- 
tive factor of order unity, but the correction to the frequency goes as rf. If the boundary conditions 
are the no-slip conditions Z — 0, then we must impose that at r = R c 

(3 — x 2 ) sinx — 3xcosa; = (A8) 
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whose first zeros are x = 5.75, 9.1, 12.3, 15.5 ... 

We have checked that the numerical results show the same qualitative dependence, even for 
polytropic or realistic stars. Knowing the dependence in the low viscosity limit allows us to regu- 
larize the equations and solve them numerically with arbitrarily small viscosity. Using s = s/rj as 
a variable, the Newtonian equations d29Tl30l become 

Z' = -Z + X (A9) 

r 

X' = --X + ( A ~^ Z- {a {a + rjs) + \ps)Z (A10) 

and can be integrated numerically since s is well behaved as rj — > 0. 

A last important remark is that for low viscosity all the overtones have very close frequencies 
(rf corrections), but the damping times scale with k 2 — (Tq w k 2 , so that modes with shorter 
wavelength are damped faster than the fundamental mode (typically a factor 4-10). Therefore, 
understanding the non-linear transfer of energy between the fundamental r— mode and higher 
order overtones may be relevant for astrophysical applications. 
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